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Abstract 


This paper introduces a general version of the information matrix 
consisting of the autocorrelation and cross- correlation matrices of the 
shifted input and output data. Based on the concept of data corre- 
lation, a new system realization algorithm was developed to create a 
model directly from input and output data. The algorithm starts by 
computing a special type of correlation matrix derived from the infor- 
mation matrix. The special correlation matrix provides information 
on the system- observability matrix and the state-vector correlation. 
A system model was then developed from the observability matrix in 
conjunction with other algebraic manipulations. This approach leads 
to several different algorithms for computing system matrices for use 
in representing the system model. The relationship of the new algo- 
rithms with other realization algorithms in the time and frequency 
domains is established with matrix factorization of the information 
matrix. Several examples are given to illustrate the validity and use- 
fulness of these new algorithms. 


1. Introduction 

Recently, system identification has gained much support for active control of flexible struc- 
tures, including acoustic noise reduction, jitter-induced vibration suppression, and fine pointing 
of spacecraft antenna. In practice, control designs based on analytical models will not work 
when used the first time because often the analytical models used in the control designs are 
not accurate enough to meet specified performance requirements. As a result, most practicing 
engineers conduct experiments to either tune the control parameters or identify accurate math- 
ematical models from input and output data. In addition to identifying system models, most 
robust control methods require information about model uncertainties. 

System identification encompasses many approaches, perspectives, and techniques (refs. 1-9) 
and can be divided into two groups. One group of techniques uses the nonparametric approach 
(refs. 1-8) with the least-squares method to determine the input-output map. The input-output 
map is characterized by a system model which may not have any explicit physical meaning. 
These techniques are generally referred to as black box approaches and are noniterative. A 
second group of techniques uses the parametric approach (ref. 9) to determine system model 
parameters. The parameters may represent physical quantities such as structural stiffness or 
mass. Nonlinear mathematical optimization techniques are used to search for the optimal value 
of each parameter. 

Both parametric and nonparametric approaches have many successful applications. While 
development and enhancement of the many individual techniques continue, a comprehensive 
yet coherent unification of the different methods is still needed. In the past decade, researchers 
have successfully provided frameworks to unify different techniques via system realization theory 
(ref. 1). 

This paper is based on the need to improve system identification techniques such as those 
discussed in references 6-8. Computational time and numerical accuracy are an issue when 
measurement data are substantial. Section 3 of this paper describes an alternate procedure used 
to perform system identification more efficiently with the concept of data correlations presented 
in references 9 and 10. A new algorithm called “system realization using information matrix” 
(SRIM) is developed. The information matrix is similar to the matrix defined in reference 9 for 



the frequency-domain analysis, but has a general form consisting of shifted input- and output- 
data correlations in the time domain or frequency domain. A special correlation matrix is 
introduced and computed from the information matrix. The special correlation matrix reduces to 
the shifted data correlation of the pulse response if the output is a free-decay response generated 
by a pulse input. To form a discrete time model, the SRIM algorithm includes several methods of 
different merit for computing system matrices, including the state matrix, the input and output 
matrices, and the direct-transmission matrix. 

The eigensystem realization algorithm with data correlation (ERA/DC) (ref. 10) uses the 
shifted data correlation of the pulse response and factors the correlation matrix via singular- 
value decomposition to obtain the system matrices. The pulse response may be obtained with 
either a pulse input or computed from input and output data with the observer/Kalman filter 
identification (OKID) technique (ref. 11). Thus, the SRIM algorithm presented in this paper 
may be considered an extension of the ERA/DC for system identification directly from input 
and output data. 

Section 4 of this paper describes how the mathematical framework developed for the SRIM 
algorithm is used to establish the relationship among different realization techniques and how 
the different algorithms may be derived from the information matrix. The basis vectors of the 
column space and the null space of the information matrix are the key elements that provide 
the link among the different realization algorithms. Matrix factorization of the information 
matrix is presented as the most important step in the theoretical development and computational 
procedure that allows unification of many different techniques. 

Each section of this paper starts with background information and ends with numerical and 
experimental examples that illustrate the validity and usefulness of the algorithms presented in 
the paper. Comparisons of algorithms are also discussed to demonstrate the merit of different 
techniques for computing system matrices. 

2. Symbols 

A state matrix, n x n 

A c continuous-time state matrix, n x n 

B input-influence matrix, n x r 

B c continuous-time input-influence matrix 

ki ith column of B, n x 1 

C output-influence matrix, m x n 

C a accelerometer output-influence matrix 

D direct-transmission matrix, m x r 

dj zth column of D, m x 1 

G(zk) frequency-response function, m x r (at frequency variable z &) 
li identity matrix of order i 

K stiffness matrix 

k time index 

kj ith spring constant 

£ data length 

M mass matrix 
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u' N (k) 

V 
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w 
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X(k) 

x(k) 

Y P (k) 


number of outputs 
ith mass 

order of first-order system 
number of zero singular values 
observability matrix, pm x n 
working matrices 

integer determining maximum order of system 

information matrix, p(m + ^) x p(m -f- r) 

fundamental SRIM correlation matrix, pm x pm 

autocorrelation of input vector, pr x pr 

cross correlation of state and input vectors, n x pr 

autocorrelation of state vector, n x n 

state-related correlation matrix, n x n 

cross correlation of output and input vectors, pm x pr 

cross correlation of output and state vectors, pm x n 

autocorrelation of output vector, pm x pm 

number of inputs 

Toeplitz matrix, pm x pr 

left singular matrix 

special force matrix at time index k,m x r 

columns of left singular matrix corresponding to nonzero singular values 
columns of left singular matrix corresponding to zero singular values 
working matrices 

output matrix containing u p (k ) up to u p (k + N — l), pr x N 

input-force vector at time index k, r x 1 

ith input force at time index k 

vector containing u(k) up to u(k -f-p — 1), rp x 1 

vector containing u(k ) up to u(k -H N — 1), pN x 1 

right singular matrix 

process-noise vector at time index k 

displacement vector 

zth element of w 

matrix containing :r(£) up to x{k N — 1), n x N 
state vector at time index k, rt x 1 

output matrix containing y p (k ) to y p {k + N — 1), pm x N 
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y{k ) output-measurement vector at time index k, m x 1 

y p (k) vector containing y(k) up to y(k -fp — 1), pm x 1 
Zfc frequency-domain variable 

a, /?, 0 parameter matrices 

ctj ith A RX-co efficient matrix, m x m (associated with output vector y(k -f i)) 

a(^) ARX-coefficient matrix, m x m (associated with output vector y(z}.) at 
frequency variable z^) 

Pi ith ARX-coefficient matrix, m x r (associated with input vector u(k + i)) 

ARX-coefficient matrix, m x r (associated with input vector u(z^) 
at frequency variable z^) 

T working matrix 

c(k) measurement noise vector at time index k 

Q ith damping term 

E damping matrix 

£ singular-value matrix 

£« diagonal matrix containing nonzero singular values 

c ith singular value 
working matrix 

0* zero-square matrix of order i 

0 i x j zero-rectangular matrix of dimension i by j 

f pseudoinverse 

Abbreviations: 

ARX autoregressive exogeneous 

ERA/DC eigensystem realization algorithm with data correlation 

FRF frequency-response function 

IDM indirect method 

OEM output-error minimization method 

OKID observer/Kalman filter identification 

SMI subspace model identification 

SRIM system realization using information matrix 

SV singular value 

3c System Realization Using Information Matrix (SRIM) 

This section presents a system realization algorithm for computing system matrices to 
characterize the map from input to output. The system matrices are the state matrix A , the 
input matrix B, the output matrix C, and the direct-transmission matrix D. The key to the 
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SRIM algorithm is computation of the information matrix, which consists of autocorrelation and 
cross-correlation matrices of the shifted input and output data. 

Section 3 starts with a description of a discrete-time state-space model and gives key 
definitions, such as the observability matrix and the Toeplitz matrix formed from the system 
matrices. The development of the state-space model realization used to compute the set of system 
matrices [.A,i3,C, D] follows. Next, the computational steps are provided for the algorithm. 
Finally, simulation and experimental examples are described. 


3.1. State-Space Model 

A deterministic, linear, time-invariant system is commonly represented by the following 
discrete-time state-space model: 


x(k 1) = Ax{k) + Bu(k) 
y(k ) = Cx{k) + Du(k) 


(i) 


where x(k) is an n x 1 state vector at time index k, u(k ) is an r x 1 input vector corresponding 
to r inputs, and y{k ) is an m x 1 output vector associated with m sensor measurements. The 
system matrices A,B,C, and D are unknown and will be determined from given input and 
output data, that is, u(k) and y{k) for Ar = 0,1,2,.. .,£. 

Equation (1) can be written for various time shifts in a matrix form as 


'ar(l) x(2) ■ 

• x(£+l)' 


A B' 

’*(0) 

*(1) . 

• x{£) 

.3/(0) 3/(1) • 

y( f -) . 


C D 

_u(0) 

“(1) • 

• u l £ ). 


If the quantities u(k),y(k) for k = Q,l,...,£, and x(k) for k — 0,1,...,^+ 1 are known, 
then A,B,C, and D may be determined from equation (2) with the least-squares technique. 
Unfortunately, the time sequence of the state vector x(k) for k — 0, 1, ...,£+ 1 is generally 
unknown and must either be estimated or measured, if possible. In practice, usually only the 
input and output sequences u(k) and y(k) for k = 1, . . ., £, respectively, are available. Therefore, 
other forms of system equations must be found for use in system identification. 


With several algebraic manipulations, equations (1) produce 


y(k) 


■ C ' 

y(k + 1) 


CA 

y(k -f 2) 

— 

CA 2 

,y(k + p - 1). 


CAP- 1 _ 



‘ D 0 


u(k ) 


CB D 


u(k + 1) 

+ - 

CAB CB D 


u(k + 2) 


_CAP~ 2 B CAP-^B CAP~ 4 B ••• D_ 


_u(k + p — 1) . 
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where p is an integer depending on the size of the system model (the dimension of A). The 
choice of p will be shown later. Let y p {k), Op, u p (k), and T p be defined as 



y(k) 


' C ' 


u{k ) 

> 


y(k + 1) 


CA 


u{k -f 1) 


y?( k ) = 

y(k + 2) 


CA 2 


u(k + 2) 



,y(k +p- 1). 


CAP- 1 _ 


* 

1 

+ 

-S£ 

J 



' D 


0 

1 




CB 

D 





T p = 

CAB 

CB 

D 





CA P ~ 2 B CA p ~ z B CA p ~ 4 B •• 

• D . 


> 


( 4 ) 


Equation (3) thus becomes 


y p {k ) — O p x(k) + T p u p (k) 


( 5 ) 


The matrix O p of dimension pm x n is commonly called the observability matrix. The matrix 
T p of dimension pm x pr is a generalized Toeplitz matrix. Note that T p is unique even though 
matrices A,B,C , and D are not unique because the system Markov parameters D and CA^B 
are unique. The use of equation (3) to develop a system model is presented in section 3.2. 

3.2. State- Space Model Realization 

The goal of state-space system identification is to determine the unknown matrices A, B, C, 
and D which are embedded in matrices O p and T p . One approach starts with computing O p 
and T p from known input and output data. With O p computed, the output matrix C is the first 
m rows of O p . Define O p (m+ 1 : pm , :) as the matrix consisting of the last (p— l)m rows (from 
(m + l)th row to the (pm)th row) and all columns of O p . Similarly, define O p ( 1 : (p — l)m, :) as 
the matrix formed with the first (p— l)m rows and all columns of O p as follows: 


Op{m + 1 : pm, :) = 


Form the following equality: 


r ca i 


‘ C ‘ 


CA 2 


CA 


CA 3 

O p (l : (p- l)m, :) = 

CA 2 


CAP - 1 _ 


CA p ~ 2 _ 

J 


( 6 ) 



' CA ' 


' C ' 


CA 2 


CA 

O p {m + 1 : pm, :) = 

CA 3 

= 

CA 2 


CAP- 1 _ 


CAP- 2 _ 


A = O p {l : (p- l)m, :)A 


( 7 ) 


The colon in place of a subscript denotes the entire corresponding row or column. The state 
matrix can then be computed by 


A = O p ( 1 : (p — l)m, :)O p {m + 1 : pm, :) 


( 8 ) 
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where f means the pseudoinverse. The integer p should be chosen so that the matrix O p (m + 1 : 
pm , :) of dimension (p — l)m x n has rank larger than or equal to n as follows: 


(p — l)m > n 


n 

p ^ 1- 1 

m 


( 9 ) 


where n is the order of the system. 


Similarly, the first m rows and the first r columns of the matrix T p (eqs. (4)) constitute the 
direct-transmission matrix D. Define T(m + 1 : (p — l)m, 1 : r) as the matrix formed by deleting 
the first m rows of the first r columns of T p as follows: 



’ CB * 


' C ' 

T(m + 1 : (p — l)m, 1 : r) = 

CAB 


CA 

CAP~ 2 B 


CAP~ 2 


B = O p (l:(p~l)m,:)B 


( 10 ) 


Equation (10) produces 


B = o\>(l : (p — l)m, :)T(m -f- 1 : (p — l)m, 1 : r) 


( 11 ) 


To determine O p and T p , first expand the vector equation (eq. (5)) to a matrix equation as 
follows: 

Y p (k) = O p X(k) + T p U p (k) (12) 

where 


X(k) = [x(&) x(£-t-l) ••• x 

[k + N-l)] 


Y p (k) = [y p (k) y p (k + 1) ••• 

y p (k + N - 1)] 



y(k) y(k + 1) 

y(k 4- — 1) 


— 

y(k + 1) y(k + 2) 

y(k 4- N) 



.y{k + p — 1) y{k+p) 

••• y(k + p + N - 2) . 


U p (k) = [u p (k) u p (k + 1) ••• 

u p (k -j- N - 1)] 



u(fc) u(k + 1) 

u(k 4- AT — 1) - 


— 

u[k + 1) u(k -f- 2) 

u(k 4- A^) 



_u(k p — 1) u(k p) 

• • • u(k 4- p-j- N — 2). 

> 


The integer N must be sufficiently large so that the rank of Y p (k) and U p (k ) is at least equal 
to the rank of O p . Equation (12) is the key equation used to solve for O p and T p and includes 
the input- and output-data information up to the data point k p N — 2. Because the data 
matrices Y p (k ) and U p (k ) are the only information given, it is necessary to focus on these two 
matrices to extract information necessary to determine the system matrices A, B,C, and D. 
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The following quantities are defined as: 


= ^Y p (k)Y^k) ’ 
Kyu = jjY p (k)u£(k) 
K m = jjU p (k)V*(k) 
n xx = ~x(k)x T {k) 

K yx = ±Y p (k)X T (k) 

n X u = ~x{k)uj{k) _ 


( 14 ) 


where N = £ — p, with £ being the data length and p the data shift. The quantities R yy , Ruu, 
and 7Z X x are symmetric matrices. The square matrices % yy ( mp x mp), 7 Z U u (rp x rp ), and 
Rxx (n x n ) are the autocorrelations of the output data y with time shifts, the input data u 
with time shifts, and the state vector x , respectively. The rectangular matrices R yu ( mp x rp), 
Ryx ( mp x n), and R X u { n x rp) represent the cross correlations of the output data y and the 
input data u , the output data y and the state vector x, and the state vector x and input data 
u, respectively. When the integer N is sufficiently large, the quantities defined in equations (14) 
approximate expected values in the statistical sense if the input and output data are stationary 
processes satisfying the ergodic property. 


When considering equations (14), postmultiplying equation (12) by U y {k) and then dividing 
the result by N will yield 

Ryu — OpRxu ~h 'YpRuu (15) 


which, ifR n }, exists, yields 


'Tp — (R y u — Op'JZ-xujR-uu 


(16) 


The matrix inverse R u ^ exists only when integers p and N are chosen so that R uv at least has 
rank rp. Similarly, postmultiplying equation (12) by Y^f (k) yields 

Ryy — OpRyx T ^ pRyu (17) 

and postmultiplying equation (12) by X T (k) gives 

Ryx — OpRxx T 'Yp'R^ lt (IS) 


Substituting equation (16) for T p into equations (17) and (18), and then substituting the 
resulting equation for R yx into equation (17) will produce 



Ryy RyuRuuRyu ~ OpRxxOp — OpRxuRuuR'xu^p 

(19) 

Now define 

Rhh — Ryy ~ RyuR U v.Ryu 

(20) 

and 

Rxx — Rxx ~ RxvR'UuR'xU 

(21) 
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Equation (19) becomes 


'R'kh — QpR'XxQj) 


( 22 ) 


Equation (22) is the key equation for determining system matrices A and C. The quantity 
'Rrjxh is determined from the output-autocorrelation matrix 7 Z yy minus the product of the cross- 
correlation matrix % yu and its transpose weighted by the inverse of the input-autocorrelation 
matrix 72. wu . The quantity 'Jthh exists only if the input-autocorrelation matrix Finn is invertible. 
The symmetric matrix 7 Z U u is invertible if the input signal u(k) for k = 1, 2, . . . ,£ is rich and 
persistent, which results in a matrix U p {k ) of full rank, that is, rp. Assume that the input 
signal u(z) for i > k is uncorrelated with the state vector x(k ) at time step k. Stated differently, 
the current and future input data are uncorrelated with the current state. In this case, the 
cross-correlation matrix 7 Z X u becomes an n x rp zero matrix and the matrix FZ X x (eq. (21)) 
reduces to Ttxx- For example, if the input u is a zero-mean, white, random Gaussian signal, 
then H'xx — FZ X x when the data length is sufficiently long (N oo in theory). 

Define 7Z as 



The matrix 72- is defined here as the information matrix and is formed by the correlation matrices 
7 Zyy, 7 Zy U , and 'Jtuu of shifted input and output data. The information matrix contains all 
information necessary to compute the system matrices A,B, C ) and D. Factoring 7 Z yields 

^ _ 'R'VV 'Ryu Ipm ^yuR'uu ^ hh Qpmxpr 171 ^pmxpr (24) 

Ifcyu Kuu] L ®prxpm Ipr J [Qprxpm FZuu j Ipr j 

where I pm (or I pr ) is an identity matrix of order pm (or pr), and 0 pmxpr ( or 0 prxpm) is a 
pm x pr (or pr x pm) zero matrix. The product of a matrix and its transpose is either a positive- 
semidefinite or a positive-definite matrix, depending on the rank of the matrix itself. Therefore, 
the matrix product on the left side of equation (24) is a positive-semidefinite or a positive- 
definite matrix. In the matrix triple product on the right side of equation (24), the left matrix 
and its transpose (i.e., the right matrix) are both of full rank. This means that 7^/* must be a 
positive-semidefinite or a positive-definite matrix as follows: 

% hh > 0 (25) 

for the case when 7 Z U u > 0 (positive definite), which is required for the existence of . 

The left side of equation (22) (the symmetric matrix 7^/*) is known from input and output 
data and the right side is formed from the product of the rectangular matrix O p of dimension 
mp x n, which is the symmetric matrix T^ xx of dimension n x n and the transpose of O p . It 
is clear that the matrix 7 Z^h must be factored into three matrices to solve for the observability 
matrix O p , and this approach is taken in section 3.2.1. 

3.2.1. Computation of A and C 

Section 3.2.1 shows two methods for computing A and C. One method decomposes the 
full matrix 7*^ and is called the full decomposition method. The other method decomposes a 
portion of Ttfrh and is called the partial decomposition method. 

3.2. 1.1. Full decomposition method. Given matrix 7 Z^h computed from the input and output 
data, the matrix decomposition method starts with factoring 72./^ into the product of three 
matrices. Singular-value decomposition is the obvious choice to perform the matrix factorization. 
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Taking the singular-value decomposition of the symmetric matrix 'Jthh yields 


K hh = UT, 2 U T = [U n 


On x rc 0 


'uT 

On^xn 0 no 


uf. 


— u Y 2 1J T 

— i ^n z -‘n u n 


( 26 ) 


The integer n 0 — pm — n is the number of dependent columns in 7 O n xn 0 is an n x n Q zero 
matrix, and 0 no is a square-zero matrix of order n 0 . The pm x n matrix U n corresponds to the 
n nonzero singular values in the diagonal matrix E n , and the pm x n 0 matrix U 0 is associated 
with the n 0 zero singular values. Combining equations (22) and (26) produces 


^hh — OpR'XxQp — Mn'&rMn 


(27) 


The last equality produces one solution each for O p and 7i xx as follows: 

Op — lifi ( 28 ) 

and 

it X x — E^ (29) 

Equation (28) implies that the pm x n matrix U n computed from the correlation matrix 7 Z^h is a 
representation of the observability matrix O p and can be used to solve for the output matrix C 
and the state matrix A with equation (8). The first m rows of U n constitute the output matrix 
C. 

Equation (29) gives the correlation % X x (eq. (21)) as the singular- value matrix E^ of the 

correlation matrix 7 For the case where the input u is a zero-mean white-noise sequence, 

the correlation / JZ XX reduces to 71 XX (eqs. (14)), which is the correlation of the state vector 
x. The diagonal nature of E^ implies that all individual elements of the state vector x are 
linearly independent and orthogonal (uncoupled). Each individual element of the state vector 
x(Ar) represents one coordinate. The importance of each coordinate can then be measured by 
the magnitude of the corresponding singular value. 

Let the diagonal matrix E n be denoted by 


E n - diagonal [<7i, <r 2 , • • •, a n ] = 


oi 

0 


0 

a 2 


L o o 


0 - 

0 

a n _ 


(30) 


with monotonically nonincreasing <jj (i = 1, 2, . . . , n), <j\ > cr 2 > • • • > a n > 0. Accordingly, 
the strength of the elements (coordinates) in the state vector :r(&) can be quantified by the 
singular values. Assume that the singular values o-j+i, • • • , cr n are relatively small and negligible 
in the sense that they contain more noise information than system information. As a result, the 
coordinates corresponding to the singular values cr 2 -_ j_i, * •• , a n are negligible compared with the 
other coordinates. The order of the system may then be reduced from n to i by deleting singular 
values cr,- + i, • • ■ , a n . 


In practice, none of the singular values will be identically zero because of system uncertainties 
and measurement noise, implying that n 0 — 0, n = pm, and U 0 — [] (empty) without singular- 
values truncation. The observability matrix Op (eq. (28)) becomes a square matrix because U n 
obtained from equation (27) is a pm x pm matrix. The equality (eq. (9)) is violated, indicating 
that equation (8) cannot be used to solve for the state matrix A. If none of the singular 
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values are zero, at least m smallest singular values must be considered as zero before using the 
full decomposition method. In other words, the last m columns of U n must be truncated and 
treated as U 0 . To overcome this problem, a different method is presented in section 3. 2. 1.2. 

3.2. 1.2. Partial decomposition method. Regardless of which integer p is chosen, the minimum 
value of n 0 must be m (the number of outputs) to make n < pm, which will then satisfy the 
equality constraint in equation (27). There is one way of avoiding any singular- values truncation. 
Instead of taking the singular- value decomposition of the pm x pm square matrix 71^, factor 
only part of the matrix as follows: 


K hh (:,l:(p-l)m)=lt??V T = {U' n U' 0 ] 


y2 

^xn 


Onxn<> 

®noXn 0 


1 

yp 

- 

yl. 


= VnKVt 


(31) 


The dimension of Tlhhi') 1 : (p — l) m ) is P m x (p ~ l)m, meaning there are more rows than 
columns. The integer n 0 indicates the number of zero singular values and also the number of 
columns of V 0 • The integer n' 0 is the number of columns of U' 0 that are orthogonal to the columns 
of U' n . For noisy data, there are no zero singular values, that is, n 0 = 0. If no singular values 
are truncated, n' 0 = m is obtained. If some small singular values are truncated, n' 0 becomes the 
sum of m and the number of truncated singular values. Stated differently, there are at least m 
columns of U' 0 that are orthogonal to the columns of U' n in equation (31). From equation (27), 
it is easy to show that 


K hh (:, 1 : (p - 1 )ra) = O p K xx O%(:, 1 : (p - l)m) = U' n l? n Vj (32) 


which yields the following equations: 


and 


o p = u' n 


(33a) 


ti xx e>T(:,l-.(p-l)m) = zlvT (33b) 

Equation (33b) does not imply that jz xx = or O p (:, 1 : (p — l)ro) = Vj if equation (33a) is 
satisfied. One important feature of this approach is that there are always enough columns of U' 0 
available for computing B and D with or without singular-values truncation. The disadvantage 
is that the singular values no longer represent the correlation of the state vector. 

3.2.2. Computation of B and D 

Similar to computing A and C, three methods are available for computing B and D. The first 
method, called the indirect method, uses the column vectors that are orthogonal to the column 
vectors of the observability matrix. The second method makes direct use of the observability 
matrix and is referred to as the direct method. The third method minimizes output error between 
the measured output and the reconstructed output. The reconstructed output is the output time 
history obtained by using the input time history to drive the identified system model represented 
by the computed matrices A, B, C, and D. 

3. 2. 2.1. Indirect method. With matrices A and C known, the input matrix B and the direct 
transmission matrix D can be computed from the Toeplitz matrix T p (eqs. (4)). To formulate an 
equation to solve for T p , the term associated with the observability matrix O p (eq. (15)) must 
be eliminated. 
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After considering equations (26) and (28), premultipiying equation (15) by UlJ and using the 
orthogonality property of U n and U 0 yield 

Uf TZyu — 14^ TplZ-uv. 


Postmultiplying the above equation by 1 Z u ^ results in 

Mj Tp — Mjltyy7l u ^ 


(34) 


Equation (34) is the fundamental equation used to solve for the input matrix B and the direct 
transmission matrix D. Note that equation (34) does not imply that T p = because Uj 

is a rectangular matrix of dimension n 0 x mp with n Q < mp. The right side of equation (34) is a 
known quantity and the left side contains the matrix T p , which is partially known to include A 
and C and partially unknown to include B and D. Therefore, the matrix T p must be partitioned 
into two parts to extract matrices B and D. 

Let T p be partitioned as 

T p - [T p (:,l : r) 7p(:,r+l:2r) ••• T p (:, (p - l)r + 1 : pr)] (35) 


From equations (4) (for definitions of T p and O p ) and equation (28) obtain 


II 

D 

M n ( 1 : (p- l)m, :)B 



Omxr 


7^(:, r + 1 : 2r) — 

D 



M n ( 1 : (p- 2)m, :)£_ 



02m xr 


T p ( :, 2r + 1 : 3r) = 

D 



Mn(l :(p-3)m, :)B _ 


(p — l)r + 1 : pr) = 

®{p — ljmxr 


(36) 


with 0 ixj being a zero matrix of dimension i x j. The product oil4^T p becomes 

Mo Tp(-, 1 : r) =Uo(:, 1 : m)D 

+ Uo(:, m + 1 : pm)U n ( 1 : (p - l)m, :)B 

UjT p (:,r+ 1 : 2 r) =Z/J(:,m+ 1 : 2 m)D 

+ U.J (:, 2m + 1 : pm)U n { 1 : (p ~ 2)m, :)B 

UjT p (:,2r+l : 3 r) =tfj(:,2 m + 1 : 3 m)D 

+ Mj (:, 3m + 1 : pm)l4 n (l : (p - 3)m, :)B 

Mo' T p (:, (p — l)r+ 1 : pr) -Uj(:,(p- 1 )m + 1 : pm)D 
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Equations (37) can be rewritten in the following matrix form: 

M 0 p ~ Mon | g (38) 

where 

U?T p (:,l:r) 

UjT p (:,r+l:2r) 

^jT p (:,2r+l :3r) 

UjTp(:,{P~ !) r + 1 : P r ) 

: m) Uj(: ) m + 1 : pm)U n {l : (p - l)m, :) 

Uj(:, m-f- 1 : 2m) Uj{', 2m+ 1 : pm)U n { 1 : ( p — 2)m, :) 

U.J (:, 2m + 1 : 3m) llj (:, 3m + 1 : pm)U n ( 1 : ( p — 3)m, :) 

Mj(:,(p-l)m-\-l : pm) O no xn 

The dimension of U 0 j is pn 0 x pr and the dimension of U 0 n is pn 0 x (m+ n). Let the right side 
of equation (34) be denoted by 

M 0 Tl — Ujn yv n u l (39) 

where Uq'ji is an n 0 x pr matrix. Equation (38) shows that U 0 j is thus given by 

M 0 n { : , 1 : r ) 

Mon( : > r + 1 : 2r) 

^o^( : ; 2r + 1 : 3r) 

Uok(->(P- X ) r + 1 : Pr ). 
and matrices B and D can be computed by 

° B =utnU oJ (41) 

The first m rows oili\ n ll 0 p form the matrix D , and the last n rows produce the matrix B. 

Equation (41) has a unique least-squares solution for matrices B and D only if the matrix 
U on has more rows than columns. Since the dimension of U on is pn 0 x (m + n), the integer 
p must be chosen so that pn 0 > (m -f n), where n 0 = pm — n, with n being the order of the 
system. For example, if p is chosen to be {p— l)m > n, then the minimum requirement for n 0 is 
n 0 = m, which indicates that the order of the system must be determined so that pn 0 > (m + n) 
is satisfied, particularly for the case where all singular values beyond a n (eq. (30)), that is, 
<r n _ . . . , crpm, are small quantities rather than zeros. 

For small n 0 , computing matrices B and D from equation (41) is efficient in time. In practice, 
the integer n 0 results from truncating small but nonzero singular values. The truncation error 
may in turn introduce considerable error in the computed results for B and D. An alternate 
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method for computing B and D without the matrix U 0 associated with the zero singular values 
is presented in section 3. 2. 2. 2. 

3. 2. 2. 2. Direct method. Instead of using U 0 to derive equation (34), the direct method 
depends on the observability matrix Op to formulate an equation to solve for matrices B and D. 
The approach used to derive the direct method is similar to the approach used for the indirect 
method. 

First, use the notation X(k) defined in equations (13) and the state equations (eqs. (1)) to 
form 


X(k 4- 1) = [*(& + 1) x(k + 2) ••• x(k + N)] 

= AX{k) + Bu' N (k) 


where u'y(fc) is defined as 


u n(^) — [“(&) u(& + 1) ••• u(k + N — 1)] 

Substituting equation (42) into equation (12) yields 


Y v {k + 1) = O p X{k + 1) + T p U p (k + 1) 

= OpAX(k) + OpBu' N (k) + TpUp{k + 1) 

u' N {k) 

Up{k + 1)_ 

= O p AX(k) + [OpB Tp]U p+ i(k + l) 


= OpAX(k) + [OpB T p ] 


From equation (12), the least-squares solution for X{k) is 

X(k) = o\[Y t {k) - T p u p (k)] 

With equation (45) for X(&), equation (44) produces 


(42) 


(43) 


(44) 


(45) 


Y p (k + 1 ) - 0 p A0\Y p (k) = —O p AO p T p U p (k) + [O p B T p ]U p+1 (k + 1 ) 

~ { [OpB Yp\ — [OpAOjiTp 0 p m x r ] } Up -i- 1 ( ^' *r 1) 

= ru p+1 (k + l) (46) 

where Qpmxr is a pm x r zero matrix and 

F — [OpB Tp] [OpAOpTp 0 pmxr] (47) 

Note that OpB is a pm x r matrix, whereas OpAOpTp is a pm x pr matrix because Tp is a 
pm x pr matrix. Since the dimensions of O p B and OpAOpTp are different, they may not be 
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directly added to become a single matrix. Similarly, T p and 0pmx7- cannot be directly subtracted. 
However, either [OpB T p ] or [O p AO p T p Opm xr] is a pm x (pr -f r) matrix. Postmultiplying 
equation (46) by Uj + i(k + 1) thus results in 


llyuirn + 1 : (p + 1 )m, :) - O p AO\H yu ( 1 : pm, :) = Y7l uu 


(48) 


where 


K V u = jjUj+i(k + l)V* +1 (k + l) 


^ W ~ 1 (^ + *) 

Now, premultiplying equations (49) by O p and postmultiplying by 71^ yields 

ol'Ryuim + 1 : (p + 1 )m, :)7t“J - : pm, i)^^ 1 = O^Y 


(49) 


(50) 


where o\,O p — I n has been applied. With input and output data given, and matrices A and 
C determined, all quantities on the left side of equation (50) are computable. The unknown 
matrices B and D are embedded in the matrix Y. Similar to equations (37) for computing the 
product oiUjTp , compute O p Y as follows: 


Olr=[B OpTp] - [AOpTp Opmxr] 


(51) 


Similar to equations (37), with the use of the alternate expression for the pm x pr matrix T p 
(eqs. (36)), one obtains 

0lY(: ) l:r) = -A0l(:,l:m)D 

+ B — AO\,{-., m+1 : pm)O p ( 1 : (p — l)m, :)B 

C9|r(:, r + 1 : 2 r) = [O p (:, 1 : m) — AO j>(:, m + 1 : 2m)]D 
+ [pj>(:, m + 1 : pm)O p (:, 1 : (p - l)m) 

— AO p (:, 2m + 1 : pm)0p(l : (p — 2)m, :)]jB 

£>J,r(:,2r+l :3r) = [C4(:,m+1 : 2m) - ^C?J(:, 2m + 1 : 3m)]D 
+ [<9p(:, 2m + 1 : pm)C? p (:, 1 : (p — 2)m) 

— AO p (:, 3m + 1 : pm)O p ( 1 : (p — 3)m, :)]j3 

C?J>r(:, pr+ 1 : (p + l)r) = Gp(:, (p— l)m+ 1 : pm)D 
Similar to equation (38), equations (52) can be rewritten in the following matrix form: 


(52) 


@pT = @pA 


D 

B 


(53) 
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where 


OpT = 


0 \ r(:, 1 : r) 
0^r(:,r+l : 2r) 
C?Jr(:,2r+l:3r) 


_0jr(:,pr+ 1 : (p + l)r)J 


-Aoj>(:,l:m) 

AOp(m -j- 1 : 2m, :)] 

[Op(:, m + 1 : 2m)— 
AC2p(2m + 1 : 3m, :)] 


[7n — 

yiC?p(:, m 4- 1 : pm)C? p (l : (p - 1 )m, :)] 

[0p( : > m 4- 1 : pm)0 p (:, 1 : (p - l)m)- 
-AC?J(:, 2m + 1 : pm)O p { 1 : (p - 2)m, :)] 

[0p(:, 2m + 1 : pm)C? p (:, 1 : (p — 2)m)- 
A0|(:, 3m + 1 : pm)0 p (l : (p — 3)m, :)] 


LoJo.tp-i) m 4- 1 : pm) 


0 


n 


Here, 7 n is an identity matrix of order n and 0 n is a zero matrix of order n. The quantity O p p 
is a pn x r matrix and 0 p j^ is a (p 4- l)n x (m -f n) matrix. Let the left side of equation (50) be 
denoted by 

OpTi = 0\ily U (m + 1 : (p 4- l)m, :)^-~J - AO\h yu { 1 : pm , :)7£“J- (54) 

Equation (50) implies that 


<V = 


0 p7e( : > 1 ; r ) 

O p n(:,r+l : 2 r) 

0^(:,2r+l:3r) 


O p Ti(‘‘,pr 4- 1 : (p+ l)r) 
Matrices B and D can then be determined from equation (53) by 


(55) 


D 
B 

The first m rows of 0^0 p p form the matrix D , and the last n rows produce the matrix B. 

Equation (56) has a unique least-squares solution for matrices B and D only if the matrix 
O p A has more rows than columns. Since the size of O p j^ is pn x (m 4- n), the integer p must 
be chosen so that pn > (m + n), where n is the order of the system. The direct method for 
computing B and D is generally more computationally intensive than the indirect method. For 
example, when n = 10 and p — 20, the number of rows in O pJ i becomes (p4- l)n = 210. However, 
it is unnecessary to use all the rows of O p j± to solve for B and D. It is sufficient to use the number 


= oJaOpT 


(56) 
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of rows in O p A and O p p so that the rank is larger than m -f n. It should also be noted that 
more rows in O pJ [ and O p p may improve the solution, particularly when considerable system 
uncertainties are present. 

The indirect and direct methods minimize the equation errors of equations (38) and (53), 
respectively. This does not imply that the output error between the real output and the 
reconstructed output is minimized. Here, the reconstructed output is the output time history 
obtained with the input time history to drive the identified system model represented by the 
computed matrices A,B,C, and D. An alternate method that minimizes the output error is 
given in section 3. 2. 2. 3. 

3.2.2. 3. Output-error minimization method. The output-error minimization method starts 
with rearranging the output equation (eq. (5)). The rearrangement is performed to formulate a 
linear equation that explicitly relates the output vector to the elements of matrices B and D. 
The least-squares solution for matrices B and D will then minimize the output error between 
the real output and the reconstructed output. 

Use equation (5) with p — N and k — 0 to obtain 

y N { 0) = 0 N x( 0) -fi T N u p ( 0) (57) 

The matrices B and D are embedded in the term Tjyu p ( 0) on the right side of equation (57). 
The method for extracting matrices B and D from equation (57) will now be shown. 

Let the column vectors in B and D be expressed as 


B = [h h hr] 1 

o = [i h ... d r ]J 


(58) 


Each column vector bj ( i = 1,2 , . . . , r) has n elements, with n being the length of the state 
vector, and each column vector {i— 1 , 2 ,..., r) has m elements, with m being the number of 
outputs. Let the vectors b and d be defined as 



The column vector 6 is the result of stacking together all the column vectors of the input matrix 
B and column vector d includes all the column vectors of the transmission matrix D. Similarly, 
let the input vector u(Ar) be explicitly written as 


u(k) = 

u 2 (k) 


. u r {k ) . 

. . . r are scalar, \ 


( 60 ) 
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With b and d (eqs. (59)) and 7]y (eq. (3)), T J yUp( 0) may now be rewritten as 



D 


u(0) 


CB D 


u(l) 

T N u p( 0) = 

CAB CB D 


«(2) 


_CA N ~ 2 B CA n -*B CA n ~ 4 B ••• D 


_u(N -1). 



Hm( 0) 1 

HmW 
Hm( 2) 

d_ -{- 

Offlxn 

CHn( 0) 

CAU{ 0) + CU n (l) 




jV— 2 ' „ 


1). 


E CA N -^U n (k) 




L k-Q -I 

where Omxrc is an m x n 

zero matrix, 

and 



— \Im u l(k) I m U2(k ) ••• I m u r (k)] | 

~ [In u l{k) In u 2(k) ''' ^n u r{k)] J 


( 61 ) 


(62) 


and I m and I n are identity matrices of order m and n, respectively. Appendix A presents a 
simple method of computing the summation term in equation (61). The matrix size of U_ m {k ) 
is m x mr and U_ n (k) is n x nr. The purpose of rewriting the expression for TjyUp(Q) is to move 
the unknown quantities B and D outside the brackets (eq. (61)). 

Substituting equations (62) into equation (57) yields 


y A r(0) = m 


where 




- c 

Hm( 0) 

Omxn 

> 





CA 

U.m( 1) 

cu n ( 0 ) 



© = 

~x(oy 

$ = 

CA 2 

Hm{ 2) 

CAWjS) + CU„{ 1) 



d 





b 


ca n ~ 1 


JV- 2 „ , „ 






Hm(N - 1) 

E CA^- h - 2 LL n (k) 
k= 0 J 

> 


(63) 


(64) 


The vector size 0 is (n + mr + nr) x 1 and the matrix size <3> is mN x (n + mr + nr). The 
unknown vector 0 can then be solved by 


0 = $ty A ,(0) 


(65) 


where f denotes the pseudoinverse. The least-squares solution 0 does not actually satisfy 
equation (63) when the system has input and output noises. However, 0 minimizes the error 
between the actual output vector y N (0) and the computed output vector y N (0) = $0. Solving 
for the least-squares solution © can be very time consuming because the number of rows in <I> 
is m times the integer N (data length). For example, the row number can be as large as 10 000 
for a system with m = 5 outputs and N = 2000 data points. 

The SRIM algorithm was developed to compute system matrices A, B , C, and D. This section 
presented two methods for computing matrices A and C and three methods for calculating 
matrices B and D. Computational steps for programming the SRIM algorithm are given in 
section 3.3. 
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3.3. Computational Steps 

To better understand the computational procedure for the SRIM algorithm, the computa- 
tional steps are summarized as follows: 

1. Choose an integer p so that p > ^ + 1, where n is the desired order of the system and m is 
the number of outputs. 

2. Compute correlation matrices 7 Zyy of dimension pm x pm, 7 Z yu of dimension pm x pr , and 
lZ U u of dimension pr xpr (eqs. (14)) with the matrices Y p (k ) of dimension pm x N and Up(k ) 
of dimension pr x N (eqs. (13)). The integer r is the number of inputs. The index k is the 
data point used as the starting point for system identification. The integer N must be chosen 
so that £ — k — p 2 > N >> min (pm,pr), where £ is the length of the data. 

3. Calculate the correlation matrix 'JZhh of dimension pm x pm (eqs. (14)), that is, = 

Hyy — 'R-yuftuu 'foyu ■ 

4. Factor TZ^ with singular-value decomposition for the full decomposition method (eq. (26)) 
or a portion of TZ^ for the partial decomposition method (eq. (31)). 

5. Determine the order n of the system by examining the singular values of 77-/^, and obtain 
U n of dimension pm x n (eq. (26)) and U 0 of dimension pm x n 0 , where n 0 — pm — n is 
the number of truncated small singular values. The integer n 0 must satisfy the condition 
pn 0 > (m + n) for the full decomposition method. For the partial decomposition method, U n 
is replaced by U’ n (eq. (31)) and the integer n 0 is the sum of m and the number of truncated 
singular values. 

6. Let U n = Op or U’ n = Op. Use equation (8) to determine the state matrix A. The output 
matrix C is the first m rows of U n . 

7. Compute U 0 ii = UjlZyyTZ^ (eq. (39)) for the indirect method, and construct U 0 n and U Q j 
(eqs. (38) and (40)). Determine the input matrix B and the direct transmission matrix D 
from equation (41) (i.e., the first m rows oiu\ n U 0 T f° rm D and the last n rows produce B ). 

For the direct method, construct O p p and O p ^ from equation (53) and solve for matrices B 

t t 

and D by computing O^Op p- The first m rows of O^O p i form matrix D , and the last n 
rows produce matrix B. 

For the output-error minimization method, construct y y(0) and from equations (64) and 
solve for matrices B and D by computing 4>tj/y(0). The first n elements of (0) form 
the initial state vector x(0), the second mr elements give the r column vectors of D , and the 
last nr elements produce the r column vectors of B. 

8. Find the eigenvalues and eigenvectors of the realized state matrix and transform the realized 
model into modal coordinates to compute system damping and frequencies. This step is 
needed only if modal parameters identification is desired. 

9. Calculate mode singular values (ref. 1) to quantify and distinguish the system and noise 
modes. This step provides a way for model reduction with modal truncation. 

The computational steps reduce to the steps for the ERA/DC method (ref. 1) when the output 
data are the pulse-response-time history. Assume that a pulse is given to excite the system at 
the time step zero. Let k — 1 in step 2. The correlation matrices 7£ yu and 7 Z U u become null, 
and 7 ^/j/j — 7Z yy is obtained. Theoretically, the formulation 'R-hh — ^yy ~ '^yuR-uu'^yu should 
not be used for computation^ of 'Rhh ^ ^nu is n °t invertible. For special cases such as free decay 
and pulse responses, 7 reduces to 7 Zyy when the integer k is chosen at the point where the 
input signal vanishes. 
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Figure 1. Mass-spring-dashpot system 


The SRIM algorithm with the computational steps incorporated with the recursive formula 
(appendix B) is more efficient computationally than subspace model identification (SMI) 
techniques (refs. 6-8). The SMI techniques require a QR (refs. 8 and 9) factorization of a 
large matrix [Up ( k ) (k)] T , followed by a singular-value decomposition and the solution of an 

overdetermined set of equations. Furthermore, the proposed method with the concept of data 
correlation permits more physical insight than the SMI techniques. 

3.4. Examples 

To illustrate the SRIM algorithm, a numerical example and an experimental example are 
given. The numerical example uses a three-degree-of-freedom mass-spring-dashpot system. 
The experimental example uses a truss structure tested at Langley Research Center. Different 
methods of computing system matrices are compared and discussed relating to system frequencies 
and dampings. 

3 . 4 - 1 - Numerical Example 

Figure 1 shows the mass-spring-dashpot system where rrij, Q, and kj (i = 1, 2, 3) are masses, 
damping coefficients, and spring constants, respectively, Wj(t) are absolute displacements for the 
respective masses, and u j are input forces. The second-order differential equation for this system 
is 


Mil) + E w + Kw — bu 


where 



mi 

0 

0 ' 


Ci + C2 

-C2 

0 

M - 

0 

m 2 

0 


-C2 

C2 + Ca 

“C .3 


_ 0 

0 

m3 


0 

-Ca 

C .3 . 


I< = 


k \ + &2 —&2 0 

— &2 &2 “b — ^3 


0 — k% ^2 + k% 

The corresponding first-order differential equation is 

x = A c x + B c u 

where 


' 

'“'1 ' 


0 

w = < 

w 2 

> b = < 

0 

. 



l 1 J 


x = 


Af- — 


O3 h 
■M~ l K -M~ l E 


B r = 


0.3x1 

6 
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with O3 being a zero matrix of order 3, I3 an identity matrix of order 3, and 0,3 X 1 a 3 x 1 
zero vector. For simplicity, let the mass mi = m2 = m3 = 1, = 1,&2 — 2, £3 = 3, and 

S — O-Ohv^- The damping matrix E = O.Obv^ is chosen so that each mode has a proportional 
damping of 0.5 percent. With a sampling rate of 1 Hz, the system is excited by a discretized 
random force u(k ) for k — 1,2,..., 3000 (normally distributed with signal of unit strength). Two 
output signals are obtained from the accelerometers located at positions 1 and 2 as follows: 



'1 

0 

o' 


*1 

y = 

0 

1 

0_ 


1 

1 


where 



0 0 
1 0 


The measurement vector y after substitution of w becomes 


y = C a M *( bu — Hii; — Kw) — Cx -f- Du 


where 

C = -C a M~ l [K E] and D = C a b 
The discrete state-space model is 


x(k + 1) = Ax(k ) + Bu{k ) + v(k) 
y(k ) = Cx(k) + Du{k ) + 6“(^) 

where 


A = 10“ 2 x 


B = 10 -2 x 


' -4.0315 

47.981 

17.092 

47.981 

-26.375 

71.972 

17.092 

71.972 

10.211 

-133.02 

19.188 

54.697 

19.188 

-70.165 

28.782 

. 54.697 

28.782 

-87.442 

r 0.72516 I 




9.6335 

39.563 

3.9627 

33.292 

.62.440 


59.137 

22.195 

3.9627 

22.195 

42.886 

33.292 

3.9627 

33.292 

62.440 

-4.8427 

47.944 

17.343 

47.944 

-26.773 

71.916 

17.343 

71.916 

9.6095 


C = 10 -2 


-300 

200 


200 0 

-500 300 


-1.6119 0.60778 0.18028 

0.60778 -1.9492 0.91167 


D = 


'0 

0 


and k is the time index. The quantities v(k) and e(k) are added to the model to represent the 
process noise and the measurement noise, respectively. In practice, both noises are generally 
assumed random-normally distributed. The process noise is set at approximately 10 percent of 
the input force, and the measurement noise is set at about 10 percent of the output force, both 
as standard deviation ratios. 

The initial values for p are arbitrarily selected as p = 6,12,25,50, and 100 to make the 
maximum system order pm = 12,24, 100, and 200, which is higher than the anticipated system 


21 



order of n — 6. For illustration, the partial decomposition method is used to compute the state 
matrix A and the output matrix C. The indirect method and the output-error minimization 
method are used to compute the input matrix B and the direct-transmission matrix D. 

In practical applications, the correct system order is unknown and the maximum system 
order selected by pm provides an upper bound, allowing the solution approach to proceed. From 
these upper bounds, there are two ways to reduce the model to the order of 6. One way is to 
truncate the singular values of 7lhh an( ^ retain only the six largest values, resulting in system 
matrices A of 6 x 6, B of 6 x 1,0 of 2 x 6, and D of 2 x 1. The eigenvalues of system matrix 
A lead to the frequencies and dampings, and the matrices C and B can be used to estimate the 
mode shapes and modal participation factors (ref. 1). 

Table 1 lists the true modal frequencies and damping ratios. Table 2 shows the identified 
damping ratios for different values of p. The identified frequencies are not shown because they are 
identical to the true frequencies up to three digits. The last column, “Error max SV,” gives the 
largest singular value of the error matrix between the real output and the output reconstructed 
from the identified system matrices with the same input signal. The error matrix has the size 
of m x £, where m is the number of outputs and £ is the length of the data. 

Another way to reduce the model order is to obtain the full-size system matrices first (i.e., 
no singular- values truncation) and then perform the modal truncation. Each identified mode is 
weighted by its mode-singular value (ref. 1, pp. 139-143). The three modes with largest mode- 
singular values are then chosen to represent the system. Table 3 shows the identified damping 
ratios for different values of p. The identified frequencies are identical to the true frequencies up 
to three digits. (See table 1.) 


Table 1. True Modal Parameters 


Mode 

True damping 
ratio, percent 

True 

frequency, Hz 

1 

0.50 

0.80 

2 

.50 

.28 

3 

.50 

.44 


Table 2. Identified Damping Ratio Obtained by Partial Decomposition 
Method With Singular Values Truncation 


P 

Mode 1 
damping, 
percent 

Mode 2 
damping, 
percent 

Mode 3 
damping, 
percent 

IDM 
error, 
max SV 

OEM 
error, 
max SV 

6 

3.06 

0.03 

1.01 

180.07 

167.5 

12 

1.13 

.53 

.50 

136.60 

120.77 

25 

.55 

.43 

.45 

128.28 

126.23 

50 

.40 

.41 

.43 

120.00 

127.04 

100 

.37 

.41 

.37 

132.00 

128.01 
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Table 3. Identified Damping Ratio Obtained by Partial Decomposition 
Method With Modal Truncation 


p 

Mode 1 
damping, 
percent 

Mode 2 
damping, 
percent 

Mode 3 
damping, 
percent 

Error, 
max SV 


3.80 

1.34 

1.12 



1.65 

.72 

.60 



.66 

.50 

.54 


50 

.46 

.47 

.50 


100 

.42 

.65 

.39 



Table 2 shows that all damping ratios are underestimated by singular-values truncation of 
'R-hh w hen the integer p is chosen sufficiently large compared with the minimum requirement 
of p , which is 3 in this example. An integer p exists that produces a minimal output error. 
Increasing the value of p beyond the optimal value does not necessarily reduce the output error. 
Both the indirect method and the output-error minimization method show similar trends on the 
output error. The output-error minimization method gives better output-error solutions than 
the indirect method, but takes three orders of magnitude longer for computation when compared 
with the indirect method. 

On the other hand, table 3 shows that some damping ratios may still be overestimated by 
modal truncation for a large p. Also, there exists an integer p that produces a minimal output 
error. Comparing tables 2 and 3 indicates that the modal-truncation method performs somewhat 
better at p = 50 than both the output-error minimization method and the indirect method. The 
computational time for the modal-truncation method is about the same as the indirect method, 
but much less than the output-error minimization method. It should be noted that the output- 
error minimization method only minimizes the output error relative to ar(0) and matrices B and 
D (given matrices A and C), and therefore, does not globally minimize the output error relative 
to all matrices A, B , C, D, and x(0). 

3. 4-2. Experimental Example 

The experimental results given in this section illustrate the usefulness of the proposed 
techniques when used in practice. Figure 2 shows the L-shaped truss structure used. This 
structure consists of nine bays on its vertical section and one bay on its horizontal section which 
extend 90 in. and 20 in., respectively. The shorter section is clamped to a steel plate that is 
rigidly attached to the wall. The square cross section is 10 in. by 10 in. Two cold air jet 
thrusters located at the beam tip serve as actuators for excitation and control. Each thruster 
has a maximum thrust of 2.2 lb. Two servoaccelerometers located at a corner of the square 
cross section provide the in-plane tip acceleration measurements. In addition, an offset weight 
of 30 lb is added to enhance dynamic coupling between the two principal axes and to lower the 
structure fundamental frequency. For identification, the truss is excited with random inputs 
to both thrusters. The input-output signals are sampled at 250 Hz and recorded for system 
identification. A data record of 2000 points is used for identification. 
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Figure 2. Truss structure test configuration. 


Table 4 lists the modal frequencies and damping ratios identified with the partial decomposi- 
tion method for determining matrices A and C in conjunction with the output-error minimiza- 
tion method for computing matrices B and D. The initial index p is arbitrarily set to make 
the maximum system order pm — 10, 20,30,40,50, 100, and 200. The singular-values truncation 
is used to reduce the order of the system model to 6. The output error decreases continuously 
as p increases from 5 to 100. The speed of decreasing the output error is slow from p = 50 to 
p — 100. The frequencies identified for all different p values are close and the damping ratios 
range from 2.5 percent to 0.4 percent for the first mode, from 1.5 percent to 0.4 percent for the 
second mode, and from 1.2 percent to 0.07 percent for the third mode. 

Results from the indirect and direct methods for computing matrices B and D are not shown 
because these methods produce output errors several orders of magnitude higher than the output 
errors shown in table 4. Both methods work well for the simulation data, with input and output 
noises assumed to be white, random, Gaussian, and zero-mean. Therefore, it is believed that 
noise nonlinearities are the major causes of the problem of introducing significant errors in 
matrices B and D. This example shows that the indirect and direct methods should not be used 
in practice for computing matrices B and D if matrices A and C are obtained by the reduced 
model with singular-values truncation. 


Table 4. Identified Modal Parameters Obtained by Partial Decomposition Method With 

Singular- Values Truncation 


p 

Mode 1 

Mode 2 

Mode 3 

OEM 
error, 
max SV 

Frequency, 

Hz 

Damping, 

percent 

Frequency, 

Hz 

Damping, 

percent 

Frequency, 

Hz 

Damping, 

percent 

5 

5.89 


7.29 

1.54 

48.5 

1.19 

621.77 

10 

5.88 


7.30 

.99 

48.2 

.71 

545.58 

15 

5.88 

1.10 

7.30 

.51 

48.0 

.99 

365.04 

20 

5.88 

.62 

7.29 

.38 

47.5 

2.06 

264.57 

25 

5.87 

.46 

7.29 

.42 

47.4 

2.44 

207.63 

50 

5.86 

.44 

7.29 

.41 

48.4 

.46 

174.25 

100 

5.85 

.43 

7.28 

.43 

48.6 

.07 

167.64 
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Table 5. Identified Modal Parameters Obtained by Partial Decomposition Method 

With Modal Truncation 


V 

Mode 1 

Mode 2 

Mode 3 

IDM 
error, 
max SV 

OEM 
error, 
max SV 

Frequency, 

Hz 

Damping, 

percent 

Frequency, 

Hz 

Damping, 

percent 

Frequency, 

Hz 

Damping, 

percent 

5 

5.89 

3.50 

7.28 

2.30 

49.0 

1.13 

817.99 

735.99 

10 

5.87 

.65 

7.29 

.47 

48.6 

.74 

262.81 

202.84 

15 

5.85 

.40 

7.28 

.41 

48.6 

.46 

197.32 

171.33 

20 

5.85 

.37 

7.28 

.41 

48.7 

.44 

216.79 

174.46 

25 

5.85 

.38 

7.28 

.42 

48.7 

.64 

203.56 

174.06 

50 

5.85 

.38 

7.28 

.44 

48.6 

.47 

198.69 

175.11 

100 

5.85 

.40 

7.28 

.45 

48.5 

.30 

194.58 

174.02 


Table 5 lists the modal frequencies and damping ratios identified with the partial decomposi- 
tion method for determining matrices A and C combined with the indirect method for computing 
matrices B and D without singular-values truncation. The full-size model is then reduced to 
the order of 6 (including only those modes of interest) and used to compute the output error. 
The output error decreases quickly when p increases from 5 to 10, and reaches a minimum at 
p = 15. The output error increases slightly again and then reduces to another minimum at 
p = 100. However, the minimum output-error value of 194.58 at p = 100 improves little from 
the minimum output-error value of 197.32 at p = 15. Similar to the results shown in table 4, the 
frequencies identified for all different p values are very close and the damping ratios range from 
3.5 percent to 0.4 percent for the first mode, from 2.3 percent to 0.45 percent for the second 
mode, and from 1.13 percent to 0.3 percent for the third mode. 

The similarity of modal parameters in tables 4 and 5 is not surprising because both tables 
share the same observability matrix before singular- values truncation. Table 4 shows the modal 
parameters computed after singular- values truncation, i.e., some columns of observability matrix 
corresponding to small singular values are truncated. On the other hand, table 5 shows the modal 
parameters computed from the full-size observability matrix. The modal parameters shown in 
table 5 are chosen to represent the system. 

The question may arise whether the output errors in table 5 may be reduced if matrices B 
and D are recalculated with the output-error minimization method with the same matrices A 
and C. The last column of table 5 provides an answer to this question. Indeed, the output errors 
are somewhat improved for all cases. The output-error value of 171.33 for p = 15 in table 5 is 
better than the output-error value of 174.25 in table 4 at p = 50. This result indicates that the 
combination of singular-values truncation combined with output-error minimization may not 
produce the global minimum for any given p value. As a result, modal truncation combined 
with the output-error minimization method works well for model reduction. 

3.5. Section 3 Summary 

A new system realization algorithm is developed with a data correlation matrix to compute an 
observability matrix with singular-value decomposition. The data correlation matrix is formed by 
the autocorrelation matrix of the shifted output data subtracted from cross correlation between 
shifted input and output data weighted by the inverse of the autocorrelation matrix of the shifted 
input data. The observability matrix is then used to compute the state matrix A and the output 
matrix C. 
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Two computational methods are presented, including a full decomposition method and a 
partial decomposition method, to determine matrices A and C. The partial decomposition 
method seems easier to use than the full decomposition matrix because the partial decomposition 
method eliminates the need for singular-values truncation. In practice, there are no zero singular 
values regardless of how clean the data sequence is. Determining the number of singular values to 
truncate requires engineering judgment or use of special techniques such as sensitivity analysis. 

Based on the computed matrices A and C, three methods are described for computing the 
input matrix B , the direct transmission matrix Z), and the initial state vector x(0). The 
indirect method uses the matrix with column vectors orthogonal to the observability matrix in 
conjunction with a data-correlation equation to extract matrices B and D. The direct method 
uses the observability matrix and the output equation to formulate an equation to solve for 
matrices B and D. The output-error minimization method determines matrices B and D and 
initial state vector x(0) by minimizing the error between the test output and the computed (i.e., 
reconstructed) output. When the input and output noises are white, Gaussian, and zero-mean, 
any combination of the methods mentioned in this summary for computing matrices A, B , C, 
D, and initial state vector x(0) performs well. For other noises, any combination works well if no 
singular- values truncation is conducted. With singular- values truncation for model reduction, 
the combination of the partial decomposition algorithm and output-error minimization works 
better than the other methods and is comparable to the modal truncation technique. 

4. Unification of SRIM and Other Methods 

Many system realization algorithms start with a state-space, discrete-time linear model and 
then formulate fundamental equations based on data correlation to compute system matrices. 
On the other hand, other algorithms use the finite-difference model and data correlation to solve 
for the system matrices. The two approaches appear fundamentally different because they use 
different types of models. The different models may yield identification results with noticeable 
discrepancy because of the different equation errors that are minimized. As a result, it is very 
difficult to interpret identification results and choose the technique best suited for a problem. 
The need exists to provide a comprehensive, yet coherent, unification of the different techniques. 

Section 4 of this paper establishes the relationship among the SRIM presented in section 3 
and other realization methods. The SRIM is derived with the state-space, discrete-time linear 
equation to form a special type of data correlation for system realization. Other realization 
algorithms, such as the observer/Kalman filter identification (OKID) technique (refs. 11 and 12), 
start by computing the coefficient matrices of the finite-difference equation, which also requires 
information from input- and output-data correlation. The similarity of using data correlation 
leads to establishing the relationship between the SRIM and the OKID method. The approach 
presented in this section provides a better way to understand and interpret other techniques, 
such as OKID and the subspace identification approach developed by other researchers (refs. 6 
and 7). 

4.1. Time-Domain ARX Model 

The finite difference model is commonly called the autoregressive exogeneous (ARX) model 
by the controls community. The ARX coefficient matrices can be computed from input and 
output data by minimizing the output-equation error (ref. 11), that is, the error between the 
actual output and the estimated output. 

The discrete-time ARX model is typically written as 

®p-iy(k + P - 1) + '& p -iy(k + p - 1) + b o fl y(Z) 

= ftp- iu(k + p — 1) + (3p~\u{k + p — 1) 4- ■ • • -f pQu(k) (66) 
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where aj for i = 1, 2, . . . ,p — 1 is an m x m matrix, (3j for i — 1, 2, . . . ,p — 1 is an m x r matrix, 
y(k ) is an m x 1 output vector at the time step k ) and u(k) is an r x 1 input vector at the time 
step k. This equation relates the output sequence y{k ) to the input sequence u{k) up to p time 
steps. 


Equation (66) produces the following matrix equality: 

<xY p (k) = (3U p (k) 


where 


a = [ofQ oci 

P = [P 0 Pi 


ctp—l ] 

Pp—l ] 


Y p (k) = 

y(k) 

y(k + 1) 

y(k + 1) 

y(k + 2) 

y(k+N-l) 
y(k + N) 


.y(k + p- 1) 

y(k +p) 

••• y(k + p + N - 2). 

U p {k) = 

u(k) 
u(k + 1) 

u(k + 1) 
u(k + 2) 

u(k -f- N — 1) 

u(k + N ) 


_u(k + p — 1) 

u(k-j-p) 

• • • u(k -f p -f AT — 2). 


(67) 


( 68 ) 


Postmultiplying equation (67) by U p {k), and noting equations (14) for definitions of 7ly U and 
7 Z uu , yield 


OcJtyu — P'R'I 


P — <YR'yu'R'uu 


(69) 


Here, the existence is assumed. Similarly, postmultiplying equation (67) by Y p (&), and 

noting equations (14) for definitions of 7Zyy and 7 Zyu, yield 

Ctllyy = mlu (TO) 

Substituting equation (69) for j3 into equation (70) thus gives 


a 



— TZyjjTZ^yTZ 


T 

yu 


= 0 


or 


allhh = 0 


(71) 


where 


^ hh ~ Hyy “ '^■yu^uu^yu 

which is identical to equation (20). Equation (71) implies that the m x pm parameter matrix a 
is in the null column space of the pm x pm matrix IZ^. In other words, any m columns of U 0 
from equation (26) or U' 0 from equation (31) may be used to construct the matrix a as follows: 

a = any m rows of Uq or 77 0 (72) 
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As a result, the m x pr matrix (3 can be solved by 


(3 - any m rows of Uj7l yn 7l u ^ or U'^7l yu 7l u ^ (73) 

In practice, no zero singular value may exist to form the matrix U 0 because of system 
uncertainties, including input and output noises. The m columns corresponding to the m smallest 
singular values may then be selected to form U 0 . 

Another way of computing a and (3 is to combine equations (69) and (70) to yield 


[-a /?] 


R-yy Ryu 

\Ryu Ruu 


— 0 


Tnxp(m-\-r) 


or 


Q3Z — 0 m xp(m+r ) 

where O^xp^.^ is an m x p(m + r) zero matrix and 


0 = [— a f3] and 71 = 


Ttyy 

7 1 

7l T 

yu 

7 1 


(74) 


(75) 


(76) 


The size of matrix 0 is mxp(m-fr) and the size of 71 is p(m+r) xp(m+r). Equation (75) implies 
that the unknown matrix 0 exists in the null column space of the matrix 7Z. Any m columns 
generated from the basis vectors of the null column space of 7Z may be considered as the solution 
for 0-^. In theory, the integer p must be chosen large enough so that the p(m + r) x p(m + r) 
matrix 7 Z is rank deficient. Because the pr xpr input correlation matrix 7 Z U u is required to have 
the full rank of pr, the rank deficiency of 7t implies that the rank of 7Z yy must be less than 
pm. Specifically, the rank of 7 Zyy cannot be more than pm — m to have at least m independent 
columns to generate the null column space of 7t of dimension m. 

4.1.1. Ob server /Kalman Filter Identification (OKID) Algorithm 

One solution for equation (75) can be derived by premultiplying equation (75) by the m x m 
matrix to obtain 

p-1 

a p— 1®7?- — 0mxp(m+r) (77) 

or, from equations (68) and (76) 

[ —Qj p-l a 0 ~~ a p- l a l a p-iP() a p ~iPl Up-iPp-l ] 77- = 0 mxp (m+r) 

Now, take out the m rows of 71 from (p — 1) * m + 1 to pm and move them to the right side of 
equation (77) as follows: 

= 7£((p — l)m + 1 : pm, :) (78) 


0 


71(1 : (p- l)m, :) 

7 Z(pm + 1 : (p + r)m, :) 


where 

0 = [-c*"Vo ~ a p-l a 1 ' • • ~ a p-l a p-2 a p- lA) a p-\Pl ■ ' • a p-lPp~l ] (79) 
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Note that the dimension of 0 is m x ((p — l)m+pr), which is m columns shorter than 0 because 
of removal of the identity matrix l m . Equation (78) should have the least-squares solution for 
0 as follows: 

© = 7£((p — l)m + 1 : pm, :) 

where f means pseudoinverse. To avoid the pseudoinverse, delete m columns of 7£((p— l)m-j- 1 : 

1(1 : (p- l)m, :) 
m + 1 : (p + r)m, :) 

The parameter matrix 0 obtained from equation (80) is identical to the parameter matrix 
used in the OKID technique (refs. 1 and 7) to identify system matrices A, B , C, D , and an 
observer gain G. The correlation matrix 7Z is referred to as the information matrix and includes 
the correlation information between input and output data. 

In theory, all methods which start with the same correlation matrix 7 Z should produce 
the same identification results. In practice, the identification results may be somewhat 
different because of the presence of system uncertainties. For example, the matrix 0 solved 
by equation (80) minimizes the output residual between the measured output and the ARX 
computed output (refs. 1 and 7). On the other hand, other methods minimize the error 
between the measured output and the output computed from an identified state-space model. 
Because different error criteria are used for minimization, identification results are expected to 
be somewhat different. Nevertheless, the results should not be significantly different unless the 
identified models are considerably reduced by either singular- values truncation and/or modal 
truncation. 

4-1.2. Experimental Example 

The data taken from the truss structure (fig. 2) are used in this section for illustration. The 
OKID method is applied to determine the system matrices A, B, C, and D. Two computational 
steps are required. The first step is to use equation (80) to compute ARX coefficient matrices to 
determine system Markov parameters (pulse response). The second step is to use the eigensystem 
realization algorithm (ERA) from the computed pulse response to realize matrices A, B, C, and 
D simultaneously. In this example, no singular-values truncation is performed in ERA. The 
ERA-identified full-size model is then reduced to the order of 6, including only the modes of 
interest. The reduced model is then used to compute the output error. 

Table 6 shows the modal frequencies and damping ratios identified with the OKID method. 
The output error decreases quickly when p increases from 5 to 10, and reaches a minimum at 
p = 15. The output error increases slightly again, and then reduces to another minimum at 
p = 50. Tables 5 and 6 have identical modal parameters. The output errors in both tables are 
very close except at p = 100, where the output-error value of 882.26 in table 6 is 4.5 times larger 
than the output-error value of 194.58 in table 5. Since the modal parameters are identical in 
both tables, the error from matrices C, B, and D may be the cause of the discrepancy. Indeed, 
when matrices B and D are recomputed by the output-error minimization method, the output 
error returns to the same level in both tables. 


pm , :) from (p — 1 )m + 1 to pm and corresponding m columns of 


7 


R(p 


7 Z(pm + 1 : (p -j- r)m, :) 


( 80 ) 
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Table 6. OKID-Identified Modal Parameters 
With Modal Truncation 


p 

Mode 1 

Mode 2 

Mode 3 

OKID 
error, 
max SV 

OEM 
error, 
max SV 

Frequency, 

Hz 

Damping, 

percent 

Frequency, 

Hz 

Damping, 

percent 

Frequency, 

Hz 

Damping, 

percent 

5 

5.89 

3.50 


2.30 

49.0 

1.13 

818.05 

735.99 

10 

5.87 

.65 


.47 

48.6 

.74 

262.81 

202.84 

15 

5.85 

.40 

7.28 

.41 

48.6 

.46 

197.22 

171.33 

20 

5.85 

.37 

7.28 

.41 

48.7 

.44 

215.97 

174.46 

25 

5.85 

.38 

7.28 

.42 

48.7 

.64 

203.21 

174.06 

50 

5.85 

.38 

7.28 

.44 

48.6 

.47 

196.70 

175.10 

100 

5.85 

.40 

7.28 

.45 

48.5 

.30 

882.26 

174.02 


4.2. Frequency-Domain ARX Model 


The frequency-domain ARX model is produced by taking the z-transform of the time-domain 
ARX model. The coefficient matrices of the frequency-domain ARX model can be obtained from 
the frequency-response data by minimizing the error between the real transfer function and 
the estimated transfer function at frequency points of interest. In theory, the ARX-coefficient 
matrices obtained from the frequency-domain approach should be identical to the matrices 
obtained from the time-domain approach. 


Let G{z k ) be the transfer-function matrix of the system described by equations (1). Consider 
the left-matrix fraction (ref. 1) 

G(z k ) - a~ l (z k )l3(z k ) (81) 


where 

a ( z k) = a o + <*\ z k l + — f &p- (p-1) 

Pi z k) = A) + Pi z k l + — Pp-i z k' [P ~ l) 


(82) 


are matrix polynomials. Every aj is an m x m real square matrix, and each fa is an m x r real 
rectangular matrix. If G(z k ) represents the frequency-response function (FRF) obtained from 
experiments, the variable z k — e j(^k/£At) ^ _ 0, 1, ...,£— 1) corresponds to the frequency 
points at Sirk/iAt, with At being the sampling time interval and £ the length of data. The 
factorization in equation (81) is not unique. For convenience and simplicity, the orders of both 
polynomials can be chosen as p — 1. 

Premultiplying equation (81) by ot(z k ) produces 


= 0{z k ) 


(83) 


which can be rearranged into 

a oG{z k ) + a\G{z k )z^ 1 + f a p _i G(z k )z^ p ^ = A) + faz k 1 + 1- P v -iz k ^ (84) 

Equation (84) is the z-transform of equation (66). With G(z k ) and z^ 1 known, equation (84) 

is a linear equation. Because G(z k ) is known at z k = e j(^k/£At) ^ _ 0, ...,£— 1), there are £ 
equations available. Stacking up the £ equations yields 

= 0 mxi (85) 
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where 



' G(* 0 ) 

G(z 0 )z^ 1 

G(zi) 
G(zi)z 1 1 

G(zi-i) ' 

G{Z£- i)^Ji 

$ = 

G( z 0 ) z q {P ~ 1) 

Ir 

G(zi)z^ ip ~^ . . . 
I r ... 

G(zi-i)z^\~ l) 

l r 


r — 1 I 

z 0 i r 

z~ l l 

z \ l t ■ • • 

Z il\ I T 


— (p—1) j 
L z 0 L r 

z~ {p ~ l h 
z \ I r • • • 

— (p— l)j 

z t-l L r J 

e = [ 

-a 0 -a i • • 

1 

R 

1 

h-* 
►— 1 

‘ * * Pp~ 1 ] 


Note that <£ is an (m -f r)p x (r£) matrix and 0 is an m x (m + r)p matrix. Equation (85) is 
a linear algebraic equation, implying that the parameter matrix 0 is in the null column space 
of <£. The null column space of $ is identical to the null column space of where * means 
complex conjugate and transpose. Postmultiplying equation (85) yields 

072- — Ojt! xp(m-f-r) (87) 

where 

77 = (88) 

The zero matrix 0 mxp(m+r) °f dimension m xp(m + r) (eq. (87)) is generally much smaller than 
0 mx £ (eq. (85)). Computing <£<$* without fully forming the matrix $ is easy because of 
special configuration. Computing 0 from equation (87) rather than from equation (85) is also 
much easier. 

Equation (87) is identical in form to equation (75) except that 72. in equation (87) is a complex 
matrix and 72- in equation (75) is a real matrix. Both equations (75) and (87) are derived to 
solve for ARX coefficient matrices a and j3. Both a and j3 are real matrices. Therefore, 72- and 
72- should have common properties. Since 72- is a Hermitian matrix, that is, 72- = 72*, the real 
part of 77 is a symmetric matrix and the imaginary part of 72- is a skew symmetric matrix. It 
becomes intuitive to suggest that the real part of 72- be considered as 72- as follows: 

72- = real part of [72-] = real part of [<£$*] (89) 

Although the matrix 72- (eq. (89)) is constructed for computing the ARX coefficient matrices, 72- 
may also be used for calculating the state matrix A and the output matrix C. Similar to the 
partition shown in equation (76), let 72- be partitioned into four parts as follows: 


77 = 


72-n 



77-12 

7^22 


(90) 


where 77n is a pm x pm square matrix, 77i2 a pm x pr rectangular matrix, and 7722 a pr x pr 
square matrix. Now consideisthe following conceptual equalities: 

77 1 1 = 77 yy 77 12 = 77 yu 7722 = 77 (91) 
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The pm x pm matrix %hh defined in equation (20) can be formed as 

— 72-11 ~ 72 - 12^22 72-12 ( 92 ) 

The matrix 7?./^ can be used with either the full decomposition method or partial decomposition 
method to determine matrices A and C. On the other hand, the matrix product can 

be used with either the direct method or indirect method to determine the input matrix B and 
the transmission matrix D. 

4-2.1. Transfer- Function Error Minimization Method 

The indirect and direct methods for computing matrices B and D do not necessarily minimize 
the error between the real transfer function and estimated transfer function when small singular- 
values truncation is involved for computing matrices A and C. Similar to the output-error 
minimization method in the time domain, the transfer-function error minimization method forms 
an equation that explicitly relates the transfer function to matrices B and D with given matrices 
A and C. 

The m x r transfer function G(z} : ) ( e T (81)) has another form of expression in terms of system 
matrices: 

G{z k ) = D + C(z k l n -A)- l B (93) 

for all zj. (k = 0, 1, . . . , i — 1). Equation (93) produces 


G(*o) I 


1m 

C^h-A)- 1 I 

G{z i) 

— 

Im 

C{z x l n - A)~ l 

-G{z e - 1). 


- Im 

C(zt_il n — A) -1 . 


(94) 


where I m and I n are identity matrices of orders m and n, respectively. Given matrices A , C, 
and matrices B and D can then be determined as follows: 




rim Cizoln-A)- 1 1 

t 

G(zq) 1 

D' 

B 

= 

Im G(z\I n — A) ^ 


G{z i) 



.Im C(z£^il n - A) 1 . 


-G(z( — 1). 


(95) 


For noisy systems, the solutions for matrices B and D do not satisfy equation (94), but 
rather minimize the error in the least-square sense between the left side and the right side 
of equation (94). In practice, only the real part of equation (94) is needed to determine matrices 
B and D. 

4-2.2. Experimental Example 

This example uses experimental data taken from a large truss-type structure. Figure 3 shows a 
Langley testbed (ref. 13) that is used to study controls-structures interaction problems (ref. 10). 
The system has eight pairs of colocated inputs and outputs for control. The inputs are air 
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8 proportional and 
bi-directional 
thrusters 

8 servo (DC) 
accelerometers 



Figure 3. NASA large space-structure testbed. 

thrusters and the outputs are accelerometers. Figure 3 depicts the locations of the input-output 
pairs. The structure was excited with random input signals to four thrusters located at positions 
1, 2, 6, and 7. The input and output signals were filtered with low-pass digital filters, with the 
range set to 78 percent of the Nyquist frequency (12.8 Hz) to concentrate the energy in the 
low-frequency range below 10 Hz. A total of 2048 data points at a sampling rate of 25.6 Hz from 
each sensor are used for identification. In this example, four FRFs from two input and output 
pairs located at positions 1 and 2 are simultaneously used to identify a state-space model of the 
system. 

The integer p = 50 is sufficient to identify as many as 50 modes (for a system of dimension 
100). A state-space model is obtained with the frequency-domain version with the system order 
truncated to 80 by singular-values truncation, and then further reduced to 78 by eliminating 
an unstable mode. Figure 4 shows the reconstructed frequency-response data (dashed lines) 
compared with the experimental data (solid lines). 

Figure 4 shows the frequency response of output 1 with respect to input 1, which represents 
a strong signal. The reconstructed FRF is obtained with the identified system matrices A, B, 
C, and D , where matrices C and D are computed by the indirect method. There are 33 modes 
(corresponding to 66 states) with damping less than 1 percent. The largest singular value of the 
transfer-function error between the test and reconstructed FRFs is 128.89. Careful examination 
of figure 4 reveals that there are noticeable discrepancies near the right end points of both the 
magnitude and phase plots. One way to fix the discrepancy is to recompute matrices B and D 
with the transfer-function error minimization method. 
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Magnitude 




Figure 4. Comparison of test and reconstructed input- 1/output-l FRF’s. Test data and input- 1 /output -1 FRF 
data reconstructed from the identified system matrices. 


Figure 5 shows the frequency response of output 1 with respect to input 1 with the newly 
computed matrices B and D. The transfer-function error is reduced from 128.89 to 65.184. 
There is clearly a trade-off, as the discrepancy between the test and reconstructed FRFs (fig. 4) 
is moved from about 12 Hz to about 8 Hz (fig. 5). 

Similar results not shown in this paper were obtained for other input/output pairs. Note that 
the frequency response of output 2 with respect to input 1 represents a weak signal. The signal 
is weak because sensor 2 is orthogonal to input 1. The results show that matching is better for 
strong signals. 
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Figure 5. Comparison of test and reconstructed input-l/output-1 FRF’s. Test data and input- 1/output-l FRF data 
reconstructed from matrices B and D. 


4.3. Section 4 Summary 

The idea of data correlation leads to relating the system realization method developed in 
this paper to other methods that are based on the discrete-time finite-difference model in the 
time domain and the frequency domain. The new state-space approach provides a way to better 
understand and interpret these methods through use of the information matrix. Indeed, the 
information matrix is the common ground for computing the system matrices. The method 
developed in this paper uses the basis vectors of the column space of the information matrix 
to compute the system matrices. On the other hand, other existing methods use the basis 
vectors of the null column space of the information matrix to compute the system matrices. In 
theory, both approaches should produce identical results. However, results may be noticeably 
different in practice because of system uncertainties and measurement noises, as evidenced by 
the experimental examples given in this paper. The mathematical unification presented in this 
paper should help users interpret the identification results obtained from different methods. 
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5. Concluding Remarks 


This paper contains three main technical contributions. First, a generalized information 
matrix is introduced consisting of shifted input- and output-data correlation matrices. Second, 
a new system realization algorithm is derived with the information matrix as the basis for 
computing system matrices. Third, several system realization algorithms are unified via the 
system information matrix. For a pulse or free-decay response, the information matrix reduces to 
the shifted output-data correlation matrix. Therefore, the classical system realization algorithms 
based on a pulse response may be considered as a special case of the algorithm introduced in this 
paper, which implies that system realization algorithm unification includes classical realization 
methods. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
December 2, 1996 


36 



Appendix A 

Efficient Computation of Output-Error Minimization Method 

N—2 

Computing Y CA N k *U_ n (k) of dimension m x nr for any integer N may be achieved by 
k = 0 

making an imaginary linear system with zero initial conditions as follows: 

xjj(N + 1 ) = Axjj(N) + ejUj(N) Xjj( 0) = 0 nx i 

Zjj(N') — Cxjj(N) i— 1 , 2 , . . . , n j — 1, 2, . . . , r 

where e* is the zth column of the identity matrix of order n, and uj(k ) is the jth element of 
the vector u(k ) defined in equation (60). At any time step N , equation (96) produces a matrix 
(such as Z{N )) of m x nr as 

Z(N) = [z n (N) 221 (A) z nl (N) ... zi r (N) z 2r (N) ••• z nr (N)] (A2) 



which gives 


N—2 



Appendix B 

Efficient Computation of Information Matrix 

Because of the nature of data shifting used to form Y p {k) and U p (k), an efficient way 
of computing correlation matrices 7 lyy, Ryu-> and 'Run exists. For simplicity without losing 
generality, consider the computation oiH yu . The product of Y p {k)U y (i k) can be written as 




' k+N-1 

E v{r)u T (r) 

r~k 

k+N 

E y( r )u T (r-l) 

r=k+l 

k+N+1 

E v{t)u t {t - 2) 

T=k~ |~2 


k+N-1 

E y{r)u T (r + 1) 
' ~ k +N 

E y(r)u r (r) 

r=k + 1 
*+;V+l 

E y(r)u T (r - 1) 

T—k + 2 


k+N-1 

E y(r)u T (r+2) 
k+N 

E y(r)u T (r+ 1) 

r=A+l 

k+N+1 

E y{r)u T (r) 

T~k + 2 


(Bl) 


From the pattern appearing in equation (Bl), the user should have no difficulty filling out all 
elements not shown. For example, all diagonal elements are identical except for their upper and 
lower limits, and the same is true for all subdiagonal elements. As a result, the second diagonal 
element may be computed from the first diagonal element by 


k+N k+N-1 

Y y(r)u T (r) = Y y{ T ) uT { T ) + y(k + N)u T {k + N) - y(k)u T (k) (B2) 

r=k+l T—k 

By induction, the diagonal m x r submatrices can be computed recursively by 
k+N+i k+N+i— 1 

Y y(r)u T (r) = Y y(r)u T (r) + y(k-{- N + i)u T (k + N + i) - y(k + i)u T (k + i) (B3) 

T—k+i + 1 T—k+i 


for i = 0, 1, . . . — 1. This recursive formula indicates that each quantity is computed from its 

previous quantity. Similarly, the first upper off-diagonal submatrix on the second m rows may 
be calculated from the first upper off-diagonal submatrix on the first m rows by 

k+N k+N — 1 

N y(r)u T (r + l)= V y(T)u T (r + 1) + y(k + N)u T (k + N + 1) - y(k)v T (k + 1) (B4) 
—k + 1 r=k 


Therefore, the first upper off-diagonal m x r submatrices can be computed recursively by 


k+N+i 

Y y( T ) uT (T+i) 

r=k+i+ 1 


k+N+i- 1 

Y. y(r)u r (r + 1) + y{k + N + i)u T {k + N + i + 1) 

T—k+i 

- y{k + i)u T {k + i + 1) (B5) 
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for i — 0, 1, . . . , p — 2. Furthermore, the second upper off-diagonal m x r submatrices are 
calculated recursively by 

k+N+i k+N+i — 1 

53 y(r)u T (r+ 2) = 53 y(r)u T (r + 2) + j/(fc + TV + i)u T (k + iV + i + 2) 

r=k+i + 1 r—k+i 

- y(k + i)u T {k + i + 2) (B6) 

for i = 0, 1, ... ,p — 3. Similarly, the third, the fourth, up to the (p — l)th upper off-diagonal 
quantities can be calculated recursively as soon as their first upper off-diagonal quantities are 
known. The first upper off-diagonal quantities are the first m rows of the product Y p (k)u2(k). 

Similarly, the lower columns may be computed as soon as the first m columns are calculated. 
For example, the first lower off-diagonal m x r submatrices may be computed recursively by 

fc+jV+i — 1 

y(r)u T (r— 1) = ^ y(r)u T (r - 1) + y(k + N + i)u T (k + N + i - 1) 

T=k-\-i+l r=k+i 

- y(k + i)u T (k + t - 1) (B7) 

for i = 1, . . . ,p— 1. The second lower off-diagonal mxr submatrices may be calculated recursively 

by 


k-\-N-\~i &+jV+z-~1 

y( r ) uT ( r - 2) = 53 y(r)u T {T - 2) + y(k + N + i)u T (k + N + i - 2) 

r=/;+5t-f-l 

- y(k + i)u T (k + i - 2) (B8) 

for i — 2, ... ,p — 1. By induction, the third, the fourth, up to the ( p — l)th lower off-diagonal 
quantities can be calculated recursively as long as their first lower off-diagonal quantities are 
determined. The first lower off-diagonal quantities in equation (B8) are the first m columns of 
the product Y p (k)u2{k). 

Computing the product Y p (k)U p (k) becomes recursive as soon as the first m rows and the first 
m columns are computed. For a long data record, that is, N » 1, this recursive procedure has 
a very efficient computing time. The computations of Y p (k)Y p (k) and U p (k)U p ( k) are similar 
to the computation of Y p (k)U^(k). Because both Y p (k)Y^ (k) and U p {k)Uj(k) are symmetric, 
only the upper half or the lower half must be computed. 
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